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We improve the steady-state ab initio laser theory (SALT) of Tiireci et al. by expressing its 
fundamental self-consistent equation in a basis set of threshold constant flux states that contains 
the exact threshold lasing mode. For cavities with non-uniform index and/or non-uniform gain, the 
new basis set allows the steady-state lasing properties to be computed with much greater efficiency. 
This formulation of the SALT can be solved in the single-pole approximation, which gives the 
intensities and thresholds, including the effects of nonlinear hole-burning interactions to all orders, 
with negligible computational effort. The approximation yields a number of analytic predictions, 
including a "gain-clamping" transition at which strong modal interactions suppress all higher modes. 
We show that the single-pole approximation agrees well with exact SALT calculations, particularly 
for high-Q cavities. Within this range of validity, it provides an extraordinarily efficient technique 
for modeling realistic and complex lasers. 
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I. INTRODUCTION 

The foundation of our understanding of lasers is semi- 
classical laser theory, in which the gain medium is treated 
quantum-mechanically and the electromagnetic fields are 
treated classically. The pioneering work of Haken [T] and 
Lamb 2J showed that the Maxwell-Bloch (MB) equa- 
tions, in which the gain medium is modeled by an ensem- 
ble of two-level atoms, successfully describe the principal 
properties of lasers, including modal thresholds, lasing 
frequencies, output power, the structure of the electro- 
magnetic fields inside and outside of the laser cavity, as 
well as dynamical effects such as relaxation oscillations 
and mode, phase and frequency locking. The only prop- 
erties that cannot be obtained from the semiclassical the- 
ory are those that depend on quantum fluctuations of the 
electromagnetic field, such as the laser linewidth, ampli- 
fied spontaneous emission and photon statistics. 

Because the MB equations are coupled nonlinear equa- 
tions in space and time, few purely analytic results could 
be obtained from the theory. Those obtained gener- 
ally relied on a number of drastic approximations: the 
mode structure was assumed to be simple (e.g. spatially 
uniform or low-order gaussian modes), the openness of 
the laser system was handled either through adding phe- 
nomenological damping to closed cavity modes or by ap- 
proximating the lasing modes as quasimodes of the pas- 
sive cavity, and the nonlinear modal interactions were 
either ignored or simplified by solving the equations near 
threshold. Where such approximations could not be 
employed, reliable theoretical results could only be ob- 
tained from brute force time-domain simulations of the 
MB equations 3 , or their multi-level generalizations. 

The last two decades have seen the emergence of novel 
laser systems based on complex resonators, driven by 
advances in microfabrication and motivated by appli- 
cations to integrated on-chip optics, as well as by ba- 
sic scientific interest. Examples are VCSELs [3], mi- 



crodisk E] , spiral [7] and (wave-chaotic) deformed disk 
lasers [SHIT] , photonic crystal lasers [HI |T3] , and random 
lasers [TU [T5] . The analytical theory existing at the time 
was not readily applicable to these complex systems; the 
random laser, in particular, poses a difficult conceptual 
challenge as the corresponding laser "cavity" has finesse 
smaller than unity, meaning that it has no isolated pas- 
sive cavity resonances. At the same time, the complexity 
of some of these structures pushed realistic simulations 
of the lasing equations to the limits of practicality. Thus, 
there was a need for a robust semiclassical laser theory 
incorporating a more accurate treatment of the cavity 
modes, including both spatial complexity and openness, 
as well as the effects of nonlinear modal interactions. 

Such a theory has been proposed by Tiireci and Stone 
[IB] . These authors employed one key approximation, 
originally introduced by Haken [TJ [TT1 [18]: the inver- 
sion is assumed to be time-independent, implying the 
absence of definite phase relationships between the las- 
ing modes (hence ignoring phenomena such as mode and 
phase locking). This stationary inversion approxima- 
tion (SIA), also called the "free-running" approximation 
[191 I20j . had been previously employed in combination 
with a third-order treatment of the non-linearity, to ar- 
rive at the Haken-Sauermann (HS) equations of multi- 
mode laser theory [TJ [2T]. About a decade ago, Man- 
del and coworkers [221 [23] also used the SIA and went 
beyond the HS equations to give an infinite-order treat- 
ment of non-linear interactions in the Fabry-Perot cavity 
with non-uniform pumping. Both the HS and Mandel ap- 
proaches neglected the effect of the openness of the cavity 
on the lasing modes (we compare our theory to the Man- 
del approach in Appendix B, having already compared 
to HS in Ref. 16 J. In contrast the current approach is 
formulated for arbitrary cavities and pump profile and 
treats the openness exactly. By seeking only the steady- 
state solutions of the MB equations, Tiireci and Stone 
|16) derived a set of self-consistent time- independent non- 
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linear equations for the lasing modes and frequencies at a 
given pump. Subsequently, Tiireci, Stone and Ge [24] and 
Ge et al. [T7] developed an iterative algorithm for solving 
these nonlinear equations, eliminated the slowly-varying 
envelope approximation, and confirmed that the resulting 
solutions agreed to high accuracy with the steady-state 
results of time-domain simulations of the MB equations. 
In 2008, Tiireci et al. [251 [55] showed that even multi- 
mode random lasing in two dimensions (2D) could be effi- 
ciently calculated using this method. We will refer to this 
approach as Steady-state Ab Initio Laser Theory (SALT) 
[2"T] . The term "ab initio" refers to the fact that the only 
inputs are the dielectric function for the passive cavity 
and a few parameters to describe the gain medium. The 
SALT method thus bridges the gap between oversimpli- 
fied analytical approaches and time-domain simulations. 
Unlike the former, it describes laser cavities of arbitrary 
complexity and openness, making no assumptions about 
the nature of the lasing modes or frequencies, or the prox- 
imity to threshold. Unlike the latter, it yields direct semi- 
analytic insights into the lasing solutions. Furthermore, 
the SALT method is in general much more computation- 
ally efficient than time-domain simulations; due to the 
elimination of the time variable, it allows for calculations 
that are impractical in the time domain due to limita- 
tions in computer speed or memory. 

In the present work, we present a significant improve- 
ment to the SALT method by introducing a new basis set 
that always contains the exact threshold lasing solution. 
The properties of the new basis functions allow us to com- 
pute the lasing solutions above threshold more efficiently 
than before. They also allow us to derive an approxima- 
tion to the full SALT for high-Q lasing cavities, which 
we abbreviate as SPA-SALT (single-pole-approximation 
SALT) which is valid well above threshold in contrast to 
the HS theory [21]. The SPA-SALT approximation yields 
solutions with negligible computational effort, once the 
threshold lasing properties are known. From this simpli- 
fied theory, we derive several analytic results for the las- 
ing behavior above threshold, including relatively simple 
formulas for the thresholds of higher lasing modes. These 
results hold to infinite order in the nonlinear modal in- 
teractions and are hence quantitatively reliable. Strik- 
ingly, these results predict a "gain-clamping" transition, 
in which higher modes are prevented from turning on 
at any pump, despite the non-uniformity of the lasing 
modes. 

The remainder of this paper is organized as follows. 
In Sec. [TTJ we review the previous formulation of the 
SALT and the solution method based on the basis set 
of constant-flux states, and describes the limitations im- 
posed by this basis set. In Sec. |III[ we present the new 
basis set and the formulation of the SALT in terms of 
this basis, and examine the efficiency of the new solution 
method. In Sec. |IV[ we derive the simplified form of the 
SALT equations arising from the "single-pole approxima- 
tion" (SPA-SALT). We then solve these equations ana- 
lytically, and demonstrate good agreement with the exact 



SALT solutions. In the appendices we derive the power 
output equations of the SALT, compare the SPA-SALT 
to the earlier Mandel approach [221 123]) and calculate 
perturbative corrections to the SPA-SALT. 

II. AB-INITIO LASER EQUATIONS 

A. The SALT Equations 

The SALT description of lasing begins with the MB 
equations for an ensemble of two level atoms interacting 
with a classical electromagnetic field: 

V 2 E + - e c (r)E+ = 4ttP+, (1) 

P+ = -i(k a - i 1± _)P+ + ^E+D, (2) 

in 

D = 7|| (Do -£>)-§: [E + (P + T - P + (E + )*] . (3) 

Here, we restrict the fields to one dimension (ID), or 
to the transverse magnetic (TM) polarization in 2D, so 
that the electric and polarization fields are scalars (the 
generalization to TE modes in 2D is straightforward). 
Their positive-frequency components are E + (r,t) and 
P + (r , t) ; in these equations, we have made use of the 
rotating- wave approximation (PWA). Note that we have 
not used the standard slowly-varying envelope approxi- 
mation, employed in most treatments to eliminate sec- 
ond time derivatives; this approximation gives no bene- 
fit in the SALT approach, and is unnecessary [T7]. We 
have taken the speed of light in vacuum c to be unity; 
wavevector and frequency will be distinguished by the 
context. D(r, t) is the population inversion, and Dq(t) is 
the pump; k a is the frequency of the gain center, j± is the 
gain width (polarization dephasing rate), 711 is the popu- 
lation relaxation rate, g is the dipole matrix element, and 
e c (r) is the cavity dielectric function, which in general is 
complex and includes the material absorption inside the 
cavity. Arbitrary cavity elements, such as mirrors can be 
represented by an appropriate choice of e c (^, although 
we will focus on dielectric cavities in our examples below. 
We assume that the E + and P + fields obey a multi-mode 
ansatz 

N 

E+(r,t) = J2*^e-* k » t , 

(4) 

P+(r,t) = Y,P^r)e-^\ 

where the indices fi = 1, 2, • • • , N label the different las- 
ing modes. The total number of modes, N, is not given, 
but increases in unit steps from zero as we increase the 
pump strength D . The values of D at which each step 
occurs are the (interacting) modal thresholds, to be de- 
termined self-consistently from the theory. The real num- 
bers k^ are the lasing frequencies of the modes, which will 
also be determined self-consistently. 
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We insert the ansatz Q into Eqs. ([l}{3]), and employ 
the stationary inversion approximation D = 0. The re- 
sult is a set of coupled nonlinear differential equations, 
which are the fundamental equations of the SALT [28 : 



D(r)=D (r) 



hp - k a + ij±J ' ^ 

N 



= o, (5) 



(6) 



$ and D are now dimensionless, measured in their nat- 
ural units e c = U^Pf^yX / (2g) and e? c = h-y^/ {Aug 2 ), and 
r„ = 7^/(7^ + — A; a ) 2 ) is the Lorentzian gain curve 
evaluated at frequency fc„. Eq. ([5| is simply a wave 
equation for the electric field mode VE^, with an effec- 
tive dielectric function consisting of both the "passive" 
contribution e c (r) and an "active" contribution from the 
gain medium. The latter is frequency-dependent, and 
has both a real part and a negative (amplifying) imagi- 
nary part. It also includes infinite-order nonlinear "hole- 
burning" modal interactions, seen in the \^ v \ 2 depen- 
dence of ([6|. In addition, we make the key requirement 
that must be purely out-going outside the cavity; it 
is this condition that makes the problem non-Hermitian. 
It is worth noting that the stationary inversion approxi- 
mation is not needed until at least two modes are above 
threshold, so ^ is exact for single-mode lasing up to 
and including the second threshold (aside from the well- 
obeyed RWA). 

Let us define a finite cavity region C, such that 



D (r) = and e c (r) = ng, r £ C. 



(7) 



Although we call C the "cavity" region, e c (r) need not 
be discontinuous at its boundary. The theory applies, 
for instance, to random lasers lacking any well-defined 
boundary [UJ HB] . For our purposes, C simply defines a 
surface of last scattering (or last amplification) , a region 
outside of which there is no dielectric nor gain material 
to affect the free propagation of waves. 
We write the external pump as 



D (f) = D F(f), reC, 



(8) 



where D is the "pump strength" and F(r) a fixed 
"pump profile" , both real quantities. The simplest case, 
F{f) — 1, corresponds to uniform pumping within the 
cavity. In general F[f) need not be uniform, e.g. if the 
pump is a finite laser spot or the gain material is dis- 
tributed unevenly. 

The lasing equation now becomes 



s (0 



l + h(r) 



= (9) 



in which h(r) = J2 U ^u\^(r)\ 2 represents the spatial hole 
burning effect. Here we have introduced the abbreviation 



Previous treatments of the SALT QH [T71 [232 US] pro- 
ceeded by inverting ^ via the Green's function to yield 
an equivalent integral equation, but for our purposes it 
is more convenient to retain the differential form. 



B. Modal Output Power 

Using Eq. ^ we can determine the unknown lasing 
frequencies fc M and mode fields ^(r). From these quanti- 
ties, all other properties associated with the semiclassical 
steady state can be derived. For instance, an important 
quantity not treated explicitly in earlier versions of the 
SALT is the time-averaged modal output power V^. This 
can be obtained in two ways. First it can be calculated 
from the asymptotic out-going fields, which are directly 
calculated in some numerical approaches [25] , Alterna- 
tively, the Poynting flux through a loop enclosing a 2D 
cavity can be converted into an area integral, which gives 
the convenient expression: 



r^D F(r) 



-Im[e(r)] |^(r0| 



7m = 7i/(fcp 



«7-0- 



(10) 



M 2tt J c " ' \ l + h(r) ""^'' J J 

A more detailed discussion and derivation of the modal 
output power is given in Appendix |A") 



Threshold Lasing Modes and Constant-Flux 
States 



The lasing equation (|9| always admits the trivial solu- 
tion ^ = 0. Below the first lasing threshold, this is the 
only self-consistent solution. As Do is gradually increased 
from zero, at some value there emerges an additional 
self-consistent solution, consisting of a single lasing mode 
^^\r). Right at threshold, this mode has infinitesimal 

amplitude, ^ (r) — > 0. Hence, the hole-burning term 
h(r) is negligible and ^ reduces to a linear equation: 

[V 2 + (e c (f) + 7^0 F(r)) kl] ¥^{r) = 0. (12) 

Note that the second term in parentheses, which we will 
refer to as e g (f), is simply the linear amplifying dielec- 
tric function of the pumped gain medium. As shown in 
the following sections, this equation has a discrete set of 
non-trivial solutions, specified by the two positive real 
numbers, {D^k^), the threshold values of the pump 
and lasing frequency. Each of these solutions would be 
a perfectly valid lasing mode at threshold for that spe- 
cific pump value, assuming that all other modes are sup- 
pressed for some reason. We refer to this set of functions 
with their corresponding frequencies as the Threshold 
Lasing Modes (TLMs). They can be thought of as the 
non-interacting lasing modes, i.e. the modes that would 
turn on in the absence of modal interactions, and their 
thresholds Dy are the non-interacting thresholds. 



4 




-12. 



12 



13 



14 



15 16 
Re[k] 



17 



19 



FIG. 1: (Color online) Trajectories of scattering matrix poles 
with increasing pump strength Do. Inset: Schematic of the 
ID slab resonator used. Its length L = 1 and index n = 1.5. 
Grey dots indicate that the pump covers the whole resonator. 
Solid curves show the pole trajectories when the gain-induced 



dielectric constant e g is given by Eq. 1 12 1, with gain parame 



ters k a = 15/L and 7x = 3/L. Different symbols lying along 
each trajectory represent different pump strengths: Do = 
(filled squares), 0.1 (open squares), 0.2 (filled circles), 0.3 
(open circles), and 0.4 (filled triangles). At D = 0, the 
poles are the resonances of the passive cavity, which all have 
the same imaginary part in this case. Stars indicate the real 
frequencies of the corresponding TLMs, which arise for 



different pump values, Do = Dq in general. In Eq. (17 1 



we associate a basis set of TCF states with each TLM; each 
TCF corresponds to adding a different gain dielectric func- 
tion to the medium, which pulls a different pole through the 
same k]P . The dashed lines show this process for these four 
poles; here we define an increasing dielectric constant e g by 
e g = sr) m , (0 < s < 1), where r\ m is the TCF eigenvalue intro- 



duced in Eq. (17 1 at the frequency of the first lasing mode. 



The dashed and solid lines for the first mode (red) coincide. 



There is another interesting interpretation of the 
TLMs. The linear wave equation (12 1 defines an elec- 



tromagnetic scattering matrix which gives the out-going 
wave amplitudes in terms of the incident wave ampli- 
tudes. The outgoing-only boundary condition implies 
that the relevant solutions correspond to poles of this 
S-matrix, i.e. eigenvectors with eigenvalue tending to in- 
finity. When D = 0, these poles are just the resonances 
of the passive cavity defined by the wave equation 



[W 2 + e c (r)k 2 ] ip(r) =0 



(13) 



with an out-going boundary condition. If the cavity is 
lossless and D$ = 0, then the corresponding S-matrix 
is unitary; otherwise it is not flux-conserving. For any 
cavity in equilibrium (i.e. lossless or absorbing) these so- 
lutions only exist for complex k, with Im(fc) < 0; hence, 
outside C, these modes grow exponentially towards in- 
finity, which means that they are not physical realizable 
[151 . When Dq > 0, the dielectric function in (12 1 is not 



merely the passive e c (^)i but includes a complex non- 
equilibrium amplifying contribution e g (f) from the gain 



medium, whose effect is to move the poles "upwards" 
towards the real axis (see Fig. [I). The non-interacting 
thresholds associated with the TLMs are the values of 
the pump that move the pole corresponding to each res- 
onance onto the real axis, making it a physically possible 
threshold lasing mode. 

In the real system, once the pump reaches the smallest 
of these thresholds the solution with Dq = D nin turns on. 
This mode then begins to contribute to the hole-burning 
term in ffify. For all higher pump values this term in- 
duces nonlinear interactions by reducing the inversion, 
raising the thresholds for the higher modes, and in gen- 
eral changing both their frequencies and spatial distribu- 
tions. Thus, above the first lasing threshold we face a set 
of coupled, nonlinear differential equations ([9]), for the 
unknown interacting lasing modes "F^ (r) and frequencies 
hp. From a practical standpoint, the most efficient way 
to solve this problem is to characterize these modes with 
a tractably small set of variables, by expanding them in 
an appropriate choice of basis functions. The original 
formulation of the SALT employed the following basis 
set: 



[V 2 + e c (r)K 2 n (k)] Vn {r, k) — 0, r € C 
[V 2 + n 2 k 2 ] <p n (r,k) = 0, fiC 



(14) 



where K n are complex and fc-dependent. The basis states 
ip n (r, k) were called the "constant flux" (CF) states, and 
they satisfy an out-going boundary condition at the cav- 
ity boundary dC [29]. Within C, they obey a wave equa- 
tion with the complex frequency K n , analogous to (13 1. 
Outside, they obey a wave equation with real frequency 
k, and are required to be outgoing at infinity. The total 
electromagnetic energy flux outside C is conserved, as it 
must be for a physical mode. 

The "constant-flux" condition outside C can be satis- 
fied by a variety of complete non-Hermitian basis sets. 
The specific CF basis Q, used in Refs. [16|l24H26]. was 



chosen because of its similarity to the equation defining 
the resonances: it differs from |l~3"| ) only by having real 
k outside the cavity. If the cavity dielectric e c is con- 
stant and the pump is uniform [F = 1), then each TLM 
is a CF state; and (ii) the complex frequency K n (k) of 
the CF state is very close to the complex frequency of a 
passive cavity resonance [THJ, [23] ■ To be precise, the CF 
frequency corresponding to a TLM is 



7.1 A) A 



*7_L 



V ■ 



(15) 



If we define K n = q n — in n (suppressing k-dependence) 
and assume that the the lasing frequency is close to the 
atomic frequency, it is easily shown [23] that 



kn 



71 + 



(16) 



which is the familiar line-pulling formula for the single 
mode lasing frequency pQ, with q^, playing the role of 
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the cavity frequency and linewidth. This emphasizes the 
relationship of the SALT to earlier theories that identified 
lasing modes with passive cavity resonances. 

When the cavity dielectric function e c and/or the 
pumping profile F is non-uniform, the TLMs are not 
given by a single CF state, and each must be written as 
a superposition of CF states. In Ref. [55], it was found 
that practical SALT calculations can be performed using 
a basis of 20 — 50 CF states. However, when the pumping 
is non-uniform, the rate of convergence of the CF basis 
set is poorer. Although the CF state definition ( 14 ) takes 



e c (r) into account, it does not include the pump profile 
F(r) as an independent parameter; effectively, these CF 
states correspond to a pump profile proportional to e c (r). 

The above drawbacks motivate us to introduce a new 
basis set for the SALT equations. These basis functions 
are still CF states, in the sense that they obey the real-fc 
out-going boundary conditions. However, their definition 
accounts for non-uniformity in both the cavity dielectric 
function and the pump profile, allowing us to assign a 
basis set to each TLM, with one of the basis functions 
exactly equal to the TLM. We will see that the nonlin- 
ear above-threshold solutions can be expanded with a 
minimum number of these basis functions, resulting in a 
marked improvement in the performance of the SALT. 

To avoid confusion, we henceforth refer to the original 
CF states (14| as Uniform Constant Flux (UCF) states, 
and the new basis states as Threshold Constant Flux 
(TCF) states. 



III. THRESHOLD CONSTANT FLUX STATES 
AND SALT EQUATIONS IN CF BASES 

A. Threshold Constant Flux States 



We define the TCF states by: 

V 2 + (e c (r) +77„(fc)F(r))fc 2 ] u n (r,k) = 0, feC 
[V 2 +n 2 fc 2 ]u n (f,A:)=0, r $. C 



(17) 



where r\ n are complex and fc-dependent, and u n (r, k) are 
outgoing with frequency k at infinity. F{r) is the spatial 
pump profile defined in (jSJ). For each k, there exists a 
discrete set {u n (r, k), T] n (k) | n = 1, 2, • • • } of solutions to 
(17 1. We refer to r\ n as the TCF eigenvalue, for reasons 



that will become clear. 

Like the UCF frequencies K n , the TCF eigenvalues 
r] n (k) are complex, not real, due to the open (non- 
Hermitian) boundary condition. One can show that 
Im[77„(fc)] < 0, which implies amplifying behavior similar 
to the condition Im[Jf„] < for the UCF states. In (17 1, 
r/ n (k)F(r) plays the role of a complex amplifying dielec- 
tric function with the same spatial profile as the pump, 
so rj n (k) physically is the scale of the amplifying dielec- 
tric constant necessary for that TCF to reach threshold 



and emit at wavevector k. As previously stated, if we 
choose k = hp , then one of the basis functions matches 



the solution \E , ^'' ) for the threshold lasing equation ([12]): 



for index n such that 



(18) 



(19) 



Note that TLMs and TCF states both satisfy linear equa- 
tions, and hence have no overall scale, so the same nor- 
malization must be assumed in (18). Thus each infi- 
nite TCF basis set is associated with one true TLM, in- 
dexed by [i. Slightly above threshold, this one TCF state 
serves as a very good approximation for the first lasing 
mode. Well above threshold, the lasing mode must be 
constructed from a superposition that includes the other 
TCF states (the lasing frequency fc M will also change 
slightly from its threshold value as the pump increases, 
and the TCF states will adjust accordingly). As noted, 
these other TCF states correspond to different values of 
e g that would also lead to lasing at k^, values that are 
not realized by the two-level gain medium of the MB 
equations. In the S-matrix picture, they correspond to 
moving a different pole through the real axis at k^i ^ , as 
indicated in Fig. [T] Higher lasing modes can likewise be 
expanded using TCF states with different k^. 

The TCF states are not power-orthogonal, but obey a 
self-orthogonality relation: 



d d r F(r) u n (f, k) u n > (r, k) = 5 n 



(20) 



We use the superscript d to indicate the dimensions of 
the system here and in the following discussion. We as- 
sume degenerate rfnS are handled, as usual, by choosing 
the basis so that (20) is satisfied. It follows that any 



sufficiently regular function having the same out-going 
boundary condition (with frequency k) can be expanded 
in the TCF basis {u n (r, k)}. For the uniform case, the 
UCF and TCF states are the same, with eigenvalues re- 
lated by 



Vn (k)=e c (K 2 Jk 2 ~l). 



(21) 



Interestingly, basis states of the UCF type were first 
defined and used by Kapur and Peierls [30] in the context 
of nuclear decay, long before their introduction to opti- 
cal physics by Tiireci et al. [16 . The /c-dependence of 
the Kapur-Peierls basis set was considered inconvenient, 
and it was largely superseded by the use of S-matrix res- 
onances, which do not form a complete basis set but are 
useful when single-pole approximations are valid (and the 
amplifying behavior at infinity is ignored) [31] ■ In our 
present situation, the appearance of internal amplifying 
eigenvalues is much more natural, for there is truly a 
gain medium within the cavity! The Kapur-Peierls (CF) 
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approach, and not the resonance approach, is thus the 
natural one for describing the laser; and with the avail- 
ability of modern computers, the fact that the basis is 
/c-dependent does not pose any serious difficulty. 



B. Threshold Lasing Conditions 



We have seen that the first TLM, having frequency 



Re[D£] (arb. units) 



„(*) 



corresponds exactly to a single TCF state 

and that the other TCF states must be in- 
cluded above threshold, even though they are not possi- 
ble TLMs for the actual system. We can find the first 
TLM by computing the TCF states and over a 

range of frequencies close to the gain center k a - For for 
a fixed choice of (rj n ,k), Eq. (19 1 will yield a complex 



(unphysical) value for Do, but when Do passes through 
the real axis at k 



, the value of r) n 



defines a TLM 



according to (18-19) (see Fig. [2]). The first lasing mode 
is then the TLM with the smallest D%. The other TCFs 

for that TLM are {u m (f, k$) \ m^n}. 

To identify which r] n will generate low threshold TLMs, 
for real D , we can rewrite (19) explicitly as 



k = K _ MVn(k)} 

Im[77„(fc)] 



D = -Im[r? n (A;)]7j_ 



with k — fcjif • From the expression in brackets in ( 23 1 , 



7-1 



(22) 
(23) 



\k — fc a | should be as small as possible — and hence, 



via (22), so should |Re(?7„)|. From the prefactor in (23), 



|Im(?7 n )| should also be small, and this condition becomes 
relatively more important than the first when 7^ is large, 
i.e. the gain curve is broad. Thus, the relevant TCF 
states are those lying within a "window" around Re(ry) « 
of width ~ 7^; within this window, states with Im(^) 
closest to zero (i.e. requiring the least gain) are favored. 
This analysis agrees with the numerical results shown in 
Fig.0 

We can also express the threshold lasing mode in terms 
of the UCF modes (fTil) 
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FIG. 2: (Color online) Relation between TCF eigenvalues rj n 
and the threshold inversion Dq . Dots show the complex TCF 
spectrum r] n (k) for a random laser. Dashed curves show the 
fc-dependence of four of these TCF eigenvalues, with arrows 
indicating increasing k. S olid curves show the corresponding 
D$(k) obtained from Eq. Q, with -y ± = fc/60. A TLM oc- 
curs when one D^(k) hits the real axis, as indicated here by 
stars. The green colored TCF eigenvalue has a small imagi- 
nary part and a real part fairly close to zero; thus it leads to 
the lowest threshold (smallest Dq), consistent with the dis- 
cussion in the text. 



for 0.25 < x < 1. Only the left half of the cavity 
(0 < x < 0.5) is pumped. The TCF state corresponding 
to the first TLM, with threshold D = 0.611, is plotted 
in Fig.|3ja), along with the UCF state making the largest 
contribution to this TLM. The TCF state is tailored to 
the pump profile and is only amplified in the pumped re- 
gion, whereas, as already noted, the UCF states have no 
knowledge of the pump profile and exhibit amplification 
within the entire cavity, including the unpumped region. 
The UCF state shown in the figure represents only 54.0% 
of the total weight in this superposition [3"2] . 

In order to reproduce the actual TLM, we must super- 
pose many UCF states to cancel the amplification in the 
unpumped region. In Fig. [3jb) we plot the lasing fre- 
quencies fc^ and (non-interacting) thresholds Dq of the 
six TLMs with the lowest thresholds obtained by solv- 



As noted in Sec. ECl for non- in § and W with 20 UCF states ' respectively. The 



uniform e c and/or F it is necessary to use a superposition 
of UCF modes: 



* M (rO = ]Ta4>„(r), 

n 



Kl - fc 2 



with k — fejL . This formulation of the SALT was used in 
Ref. [35] to analyze 2D random lasers. 

To illustrate the advantage of the TCF basis for non- 
uniform cavity dielectric function e c and pumping profile 
F, we study a ID resonator of length L = 1. The re- 
fractive index is n = 1.5 for < x < 0.25, and n = 3 



£ Jj d rF{f) Vn {r)^ n ,{^<, 



(24) 



largest deviation between the TCF and UCF thresholds 
is 0.68%, and the frequency differences are below 0.1%. 

In more complex lasers, e.g. the 2D random lasers of 
Ref. [3S] , a still larger UCF basis set is required to achieve 
results comparable with the TCF basis. In Fig. |2J the 
cavity C is defined by a disk of radius R = 1 , in which we 
randomly place 600 dielectric particles of radius ~ i?/80 
and index n — 1.2. In Fig. pta), we subject the entire 
cavity to a white noise pump, 



F(r) = l+£(r), 



(25) 



with max|f(r)| = 0.3. We find that 50 UCF states 
must be included in the UCF expansion in order to 
achieve good agreement between the threshold solutions 
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(a) 



x1(f 



14 U 
Threshold lasing frequency 

FIG. 3: (Color online) (a) Spatial profile of the first thresh- 
old lasing mode of a ID slab resonator of of length L — 1 
(solid blue curve). This TLM corresponds exactly to a TCF 
state. The matching UCF state, having the largest overlap 
with this lasing mode, is shown for comparison (red dashed 
curve). Inset: Schematic of the resonator. The refractive in- 
dex is n = 1.5 for < x < 0.25, and n — 3 for 0.25 < x < 1. 
The gain center is k a — 15/ L, and only the left half of the 
cavity is pumped (indicated by the grey dots). Both states 
are normalized to unity at x = 0. (b) Lasing frequencies hp 
and non-interacting threshold values Dg of the six TLMs with 
the lowest thresholds. Crosses are given the TCF solutions, 
and squares are the UCF solutions with 20 UCF states. 



of |24| and the TCF predictions (|22| and {23}. When 
we pump only part of C (keeping the scatterer configu- 
ration fixed), more UCF states are needed to correctly 
reproduce the TLMs, even in the absence of the pump 
noise. In Fig. gb), the pump covers a central area of 
radius R/2. We find that a superposition of 200 UCF 
states is required to generate a TLM whose false-color 
intensity plot (inset) is indistinguishable by eye from the 
corresponding TCF state (not shown). Even with this 
many UCF states, the computed TLM intensity pro- 
file still differs significantly from the exact (TCF) pro- 
file when plotted along any arbitrary direction, as shown 
in the main figure. The reason so many UCF states are 
required is that the TLM (and the corresponding TCF 
state) is amplified only up to the boundary of the pump 
region; whereas each UCF state, like the uniform pumped 
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FIG. 4: (Color online) Comparison of TCF and UCF results 
for 2D random lasers, (a) Lasing frequencies and non- 
interacting threshold values Dq of the twelve lowest TLMs, 
for a white-noise pump covering the whole cavity. The pump 
profile is given by Eq. (25 I. Crosses show the TCF solutions 
for (22 1 and ( |23[ ), and squares show the UCF solutions for 
(241 with 50 UCF states, (b) Spatial intensity profiles for 
the first threshold lasing mode of a partially pumped random 
laser, with F = 1 for r < R/2 and F = for r > R/2. The 
intensity is plotted along the line 9 = 225°. The solid curve 
shows the TCF solution, and the dashed curve shows the UCF 
solution computed from 200 UCF states. Inset: false-color 
intensity plot from the superposition of 200 UCF states; the 
TCF intensity plot, which is not shown, looks similar. 



system, is amplified up to the boundary of C [35]. Us- 
ing these 200 UCF states, the calculated mode threshold 
and frequency are D = 0.142 and k = 30.011; the exact 
TCF results, from ([22} and ((23), are D = 0.140 and 
k = 30.006. 



C. Above-threshold Lasing Modes 

Above threshold, each lasing mode can be efficiently 
expanded as a superposition of TCF states: 



(26) 
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This expansion automatically satisfies the appropriate 
out-going free wave equation outside C. By inserting 
the above expansion into we write the latter as 



1 + h(r) 



7m 



(27) 



Following the procedures used in Ref. [28], we multiply 
both sides of (27 1 by F(r)u n '(r), integrate r over C, and 



invoke the self-orthogonality property (20), to find the 
SALT equation in the TCF basis: 



7 M f d< i r F(r)u n (r)u n -(r) 



Vn JC 



(28) 



l + h(r) 



Eq. ( 28 ) is a set of nonlinear fixed-point equations above 



threshold, one for each lasing mode. In general, the com- 
plex matrix 7~ nn i (k), which we refer to as the lasing map, 
has complex eigenvalues. Because the pump strength Dq 
is a real variable, the unknown lasing frequency must 
be such that one of its (nonlinear) eigenvalues is real 
and equal to I/Dq. This is achieved by tuning fc M to 
find the values at which the different eigenvalues corre- 
sponding to the different modes cross the real axis, as 
follows. (This procedure is the same as for the UCF ba- 
sis, and was described in Ref. [2S]). The first threshold 
and lasing frequency are found simply by solving (22) 



and (23) self-consistently, as described in the previous 



section; these equations are the diagonal form of ( 28 ) at 
threshold. The associated TLM is proportional to this 
solution, with vanishing overall amplitude. We then in- 
crease Dq in small increments, and use the solution for 
the smaller pump value as a starting point for the non- 
linear solver. At each step, the nonlinear solver adjusts 
the frequency k^ so that the corresponding eigenvalue of 
Tnn'ikfj,) is real. From the modified lasing map, which 
includes the hole-burning term, we can determine if a 
second mode has reached its (interacting) threshold [26] . 
A similar procedure works for third and higher modes, 
and has been shown to work for systems as complex as a 
2D random laser with eight modes turned on [25] . 

In earlier works, the lasing map was written in the UCF 
basis. This has the same form as (28), with a slightly 
different matrix operator: 



7m K 



K 2 



h.2 
M 



^ f(r >" ( :>"' (ri (29) 

1 + h(r) 



At threshold (h — > 0), we recover the threshold lasing 
equation ( 24 ) . The solution algorithm is identical to that 



for the TCF map, except that the full matrix solution 
must be performed even at the first threshold since the 
UCF map is not diagonal. 



Fi g. |5| compares the lasing modes obtained from ( 28 ) 
and ( |29| f or the ID slab resonator that we studied earlier 
in Fig. 13] For Dq = 1.264, there are two lasing modes. 



(This value of Dq is approximately twice the first lasing 
threshold, Dq~ = 0.611.) Using 20 basis functions for 
both methods, we find good agreement in the predicted 
spatial profiles. Fig.jf^b) shows the largest expansion co- 
efficients of the two modes in the TCF and UCF bases. 
We find that both modes retain a dominant component 
in the TCF basis, even when the system is significantly 
above threshold. As the pump strength increases, the 
spatial hole burning term changes e g (r), so the weights 
of the dominant components in the TCF basis gradually 
decrease, but they remain larger than 80% in the calcu- 
lated range. In contrast, the largest components of the 
two modes in the UCF basis are less than 60%. 

We remark that we could in principle absorb the hole- 
burning denominator l/[l + h(r)], calculated at the pump 
strength (Dq — 6Dq), into the profile function F(r), to 
produce an even better set of modified TCF states for 
the nonlinear calculation at Dq. This is essentially an 
alternative means of solving the non-linear problem by 
keeping the self-consistent equation almost diagonal in an 
evolving basis; however, it is usually too computationally 
expensive to recompute the TCF states this way. 



IV. APPROXIMATE ANALYTIC SOLUTION OF 
THE SALT EQUATIONS 

A. Alternative Fixed-point Equation 

Analyses of the MB equations, either in the single- 
mode or multi-mode lasing regime, almost always employ 
the near threshold approximation, in which the infinite- 
order nonlinearity of Eqs. (5| |6|28 ) is truncated at cubic 
order to give a near-threshold approximation to the so- 
lution. (An exception to this is the work of Mandcl and 
coworkers [221 [S3] discussed in Appendix B). Based on 
this cubic approximation, and the approximation of a 
closed cavity, Haken and Sauermann (HS) long ago de- 
rived a set of constrained linear equations for the modal 
intensities in the multimode regime [21] , The HS equa- 
tions have been studied further [TJ [33J, and have been 
used to analyze random and complex lasers in recent 
years [34l [35] . However, the results are unsatisfactory, 
as shown by Tiireci et al. [16] . The cubic nonlinearity in 
the HS equations leads to a saturation of modal inten- 
sities, in disagreement with the linear increase expected 
on general grounds, and found by more exact treatments 
[T71 HM] . It also allows many more modes to turn on than 
in the more exact treatments [T71 [35] • These failures 
are unsurprising, as the HS theory is being applied to 
a regime well above the first threshold, where the cubic 
approximation is poor. Generalizing the equations from 
cubic to higher orders rapidly becomes unmanageable, 
since higher powers of the intensity generate many more 
interaction coefficients to take into account. 

In Ref. [16] , it was shown that in the limit of large hole- 
burning, h(r) 3> 1, the SALT correctly predicts mode in- 
tensities growing linearly with Dq, within the single-pole 
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FIG. 5: (Color online) (a) Spatial profiles of two lasing modes 
above threshold in a ID slab resonator. Inset: Schematic of 
the resonator. The properties of the resonator are given in 
the caption of Fig. [3] The pump strength Do = 1.264 is 
slightly higher than twice th e fir st threshold. The solid lines 
and circles are the results of (281 and (129]), respectively, both 



using 20 basis functions, (b) Weights of the largest expansion 
coefficients of both modes in the TCF (solid curves) and UCF 
(dashed curves) bases. The second mode has an interacting 
threshold D = 0.89. 



approximation to be discussed below. In this current sec- 
tion, we will derive an alternative lasing map which, in 
the same approximation, allows a more complete analytic 
solution that demonstrates linear behavior for all values 
of h(f) and for multi-mode lasing. This approximation 
also provides a quantitative solution for the modal inten- 
sities, slopes and interacting thresholds, in good agree- 
ment with the exact solutions of the SALT equations. 
To our knowledge, these are the first results of this type, 
valid for arbitrarily complex cavities, to appear in the 
literature. 

In order to develop the desired approximation, we first 
re-express the lasing equations in terms of the inverse of 
the map T nn ' (k) defined in ( p8| . This inverse map has the 
same fixed points, but is much more convenient to work 
with. Wc multiply both sides of (27 1 by [1 + h(r)], and 
repeat the steps leading to (28 1, i.e. projecting the two 



7m 



(30) 



(31) 



where h nn > (k) = J c <fV F{r, k) h(r) u n (r, k) u n > (r, k) 



Note that ( 30 ) has the same form as ( 28 1 , but with the 



quantity Do, which plays the role of the eigenvalue, in- 
verted. This implies that r = T~ , which can be con- 
firmed by multiplying the two operators and using the 
completeness and self-orthogonality of the TCF states. 

The operator r, through the term h nn i(k), contains 
only a second-order dependence on the lasing modes, in 
contrast to the infinite-order dependence occurring in T . 
Thus, (30) possesses only a cubic nonlinearity, but this is 



not the same cubic nonlinearity that appears in the HS 
theory. No Taylor expansion has been performed; the 
inverse lasing map is exact at all pump values, and we 
are still working with infinite-order nonlinearity in the 
conventional sense of using a dielectric function which 
contains the field to infinite order. 

We could, in principle, use the inverse map r to solve 
the SALT equations, in the same way that we used T . 
Preliminary investigations show that such an approach 
is possible, and may have some interest, but we will not 
pursue this further here. Our aim is instead to introduce 
the "single-pole approximation" (SPA) into (30). This 



sides onto the TCF basis and using the self-orthogonality 



gives a simple approximate solution that is very easy to 
implement, and yields important analytic results. 



B. Single-pole SALT Equations 

The single-pole approximation was introduced in [16] 
to show the connection between the SALT equations, 
which solve the MB equations with minimal approxima- 
tions (principally the RWA and the stationary inversion 
approximation) |17j . and the HS equations which employ 
many more approximations. Aside from the aforemen- 
tioned cubic approximation, the HS theory assumes that 
the lasing mode is accurately described by a passive cav- 
ity mode. As we have seen, even the threshold lasing 
mode is not a passive cavity mode: it is neither a closed 
cavity mode (as assumed by HS), nor a passive cavity 
resonance as often assumed in the literature. Further- 
more, above threshold the nonlinearity mixes in other 
TCF states, which changes the spatial distribution, am- 
plitude, and frequency of the lasing mode. This effect 
is quite important in low-Q cavities, such as the ran- 
dom lasers treated in Ref. [25] , and the full SALT theory 
describes this effect very well. In high-Q cavities, the 
mixing in of other TCFs is much weaker, because the 
scattering from the gain medium is so much weaker than 
the scattering from the cavity itself. Therefore, it is rea- 
sonable to assume that the lasing modes above threshold 
have the same spatial profile as the TLM, with an am- 
plitude that can increase with Dq. This is equivalent 
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to taking only one term in the expansion of the cavity 
Green's function in the CF basis; since each term has a 
single pole in the complex plane, Tiireci et al. [16J called 
this approach the single-pole approximation (SPA). 
To be precise, the SPA assumes that 



a'Zu, 



10 (7; ) 



(32) 

where u no (r) is the TCF which is equal to the TLM at 

threshold and Kjy is the threshold value of the lasing 
frequency. With this approximation the additional index 
n is redundant and can be omitted, as we do henceforth. 
Thus the SPA assumes both that the lasing modes are 
fixed as TLMs, and that the lasing frequencies are fixed 
to be their threshold values. The remaining quantities 
to be calculated are just the number of modes and their 
amplitudes a M (Do) at a given pump value [55] . This also 
necessitates finding the interacting thresholds Dq int . 
With this approximation the nonlinear matrix equa- 



tion (301, after canceling a common factor a^, is linear 



for the modal intensity 1^ 



Do 



- 1 = 



Xflb 



^^vX^vh (33) 



d d r F(r) u 2 Jr) \u v {r)\ 2 . (34) 



Here Dq = 17^/7^ are the non- interacting thresholds 
for the TLMs, which are obtained together with fc/f at 



threshold using (22) and (23). Because the frequencies 



of the modes are assumed to be fixed, the spectral gain 
factor r„ and the "interaction constants" x^v are pump- 
independent quantities. Note also that every quantity in 
Eq. (33) is real except for Xnv, wn ich must have some 
imaginary part if the cavity is open. This inconsistency 
is a consequence of the single-pole approximation; how- 
ever, the higher Q the cavity, the smaller is the imaginary 
part, and for most cavities of interest it is acceptable to 
neglect this imaginary part. Henceforth we will use the 
approximation \nv ~ R- e [X/ii/]i anc ^ simply denote the 
real part with the same symbol. With this approxima- 
tion, the matrix Xnv is real and has positive elements. 

The above result, which we will term the SPA-SALT, 
bears a remarkable resemblance to the HS equations. 
Those equations take the form 



1 



Do 



E 



(35) 



It can be shown that the cavity decay rate k^, a quan- 
tity inserted by hand in that theory, is simply Dq in the 
SALT, which is calculable once the cavity and pump pro- 
file are given. The coupling matrix in the HS equations 
has exactly the same form as in the SPA-SALT, except 
that HS used closed cavity modes (not the real part of 
of the open cavity TLMs) , and did not take into account 
the pump profile, F(r). However, the different depen- 
dence of ( 35 1 on Dq in comparison to ( 33 ) leads to very 



different behavior at large pump values. For pumps near 
the first threshold, the two equations are approximately 
the same, but at large pump it is easy to show that the 
modal intensities in the HS theory saturate to a constant, 
whereas in the SPA-SALT they are proportional to Do. 
(It should be noted that Do in the MB equations, which 
we refer to as the pump, is actually the equilibrium value 
of the inversion in the absence of laser emission. When 
one has a multi-level laser with a true pump between up- 
per and ground levels which are distinct from the lasing 
transition, the quantity Do is a function of the pump 
which is linear at small pumps, but saturates eventually, 
and is bounded by the value corresponding to complete 
steady-state inversion of the lasing levels.) 



C. General Solution of the SPA-SALT Equations 



Let us rewrite the SPA-SALT equation (33) as 

TyT — 1 = Alii/ Ivi = Y v X[iv- (36) 

This seems to be simply an inhomogeneous linear sys- 
tem to be solved by inversion, but in fact it is more 
complicated, for we have not indicated the number of 
modes to be summed over. Let us suppose that we have 



solved the non-interacting threshold conditions ( 22 ) and 



( p3| for a given e c (r) and F(r), obtaining a subset of M 
TLMs {u^(r) | /i = 1, 2, • • • M}, with real non-interacting 
thresholds Dq less than some cut-off value, Do, c (taken 
to be much higher than the first lasing threshold). For 
a given Dq, the indices //, v occurring in (36) are those 



corresponding to lasing modes that have turned on. We 
have used this fact in deriving ((Ml), where we divided 



out the common factor a„, which is valid only if a„ is 



non-zero. Hence (33) is a constrained inversion problem; 



for each value of D , we must construct the matrix 
from the correct subset of the M TLMs at our disposal. 

We wish to find an ordered set of matrices 
A$,A$, . . . A^ ma *-\ as well as the associated interact- 
ing thresholds Dq • t , which are the values of Do at which 
the /x-th mode turns on. Because the SPA-SALT includes 
the effects of nonlinear modal interactions, these differ 
from the non-interacting thresholds Dq. In fact, iV max 
often is less than M, since some of the candidate modes 
may never turn on at any pump value, as we will see be- 
low. For a given Do, let us suppose that N lasing modes 
have turned on. Without loss of generality, we assume 
that the indices for these lasing modes are /x = 1, • • • ,N. 
We now have a non-sparse N x N matrix A^ Vl and can 



invert ( 33 1 to obtain 



I it = c, t D - b M 



H = 1,2,...N 

N 



(37) 



From this, we see that the intensity of each lasing mode 
increases linearly with Dq, between each threshold, no 
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matter how many modes are lasing or how far the laser 
is above threshold. 



(N+l) 

A 1 for the remaining set of 



To find the next matrix A^l^ 1 ' we must find the low- 



est interacting threshold D Q inj 
M — N modes. To do this, we note that (37) is valid 



for Dq- itA < D < D ^l (the lasing intensities are con- 
tinuous at each threshold although their slopes are not). 
At the upper limit of this range, Dq — D Q ^, we can 
equally well add mode (N + 1) to this matrix equation. 
The resulting equation would yield identical solutions for 
,In-i phis the solution In+i = 0. Thus we can 
evaluate ( 36 ) for all choices fi = N + 1, . . . M: 



^O.int — U Q 



N 



(38) 



which gives N—M explicit linear relations for the possible 
N + 1 st threshold. Evaluating these relations, one simply 
chooses the lowest value, which is then the correct N+l st 
interacting threshold. This defines a recursive procedure 
to find all the interacting thresholds and uniquely deter- 
mine the ordered set of A matrices required to compute 
I/j.(Do) for the entire desired range of D Q . 

Note that we always assume the "non-trivial zero" so- 
lution at each (interacting) threshold, i.e. that the physi- 
cal solution switches from the trivial zero for In+i to the 
non-zero lasing solution, giving rise to a bifurcation with 
discontinuous slope. When this happens, all the modes 
which are already turned on experience a negative kink in 
their slopes at higher thresholds. This behavior is char- 
acteristic of lasers when higher modes turn on, and the 
SPA-SALT captures it in a simple manner. 



Once the constraints on Eq. ( 36 ) are implemented in 



this manner, the solution of the SPA-SALT equations re- 
quires just N max inversions of relatively small matrices 
generated from the input parameters, {x^}, {r„}, {-Dq}. 
Thus the computational time for solving the SPA-SALT 
equations is negligible once the TLMs have been cal- 
culated. When the single-pole approximation is good, 
the nonlinear multimode problem becomes only mini- 
mally harder than the linear TLM problem, which can 
be adapted for efficient solution using finite element or 
boundary element methods [33 [35]. In Section IV F 



we compare the SPA-SALT lasing solutions to the ex- 
act SALT calculations, finding good agreement. Note 
that it has already been shown [T7] that the exact SALT 
solutions agree to within a few percent with exact time- 
dependent MB simulations for simple ID edge-emitting 
lasers, as long as the conditions for the stationary inver- 
sion approximation are well-satisfied. 



D. Gain-clamping Transition 



Eq. (37) gives a linear relation determining each of the 
N m ax interacting thresholds of the form 

1 



1-A, 



D%, (39) 



where the function / M = (1 — A M ) _1 , is the threshold en- 
hancement factor which increases the fi th threshold from 
its non-interacting value, due to the spatial hole-burning 
of lower threshold modes, which depletes the gain. In 
simplified treatments of the laser rate equations, in which 
the cavity mode is assumed perfectly uniform in space, 
these interactions actually clamp the effective gain so 
that it no longer increases with the external pump, pre- 
dicting that no additional modes turn on jj). In reality 
all resonators admit multiple modes with incomplete spa- 
tial overlap and so this extreme gain clamping behavior 
is not realized. The SPA-SALT gives a much more rig- 
orous criterion for gain-clamping at the level of the iV th 
lasing mode. If Ajy — > 1 then all higher thresholds are 
pushed off to infinity and no more modes can turn on for 
any value of the pump. 

Note the analogy here to mean-field phase transitions, 
for example where a strong enough magnetic interaction 
causes the susceptibility to diverge. Here strong interac- 
tions, meaning large values of the coefficients x^v{^ 7^ v), 
suppress "ordering" of higher modes. Conversely, spa- 
tially disjoint or weakly overlapping modes will not be 
suppressed and their interacting threshold will be ap- 
proximately equal to their non-interacting thresholds. In 
addition, higher modes with substantially lower modal 
gain and Q-values with respect to the first mode(s) will 
be more easily suppressed. Calculations for various ex- 
amples indicate that this gain-clamping "phase" of the 
laser can be reached for realistic lasers. We calculate and 
discuss the coefficient A 2 below. 



E. One- and Two-mode Solutions 

To get a feeling for the SPA-SALT solutions, we now 
present explicit results for one and two mode lasing, 
which illustrate most of the qualitative features of the 
theory. The single mode result is trivial. The lowest 
non-interacting threshold, Dq , is found as part of the 
calculation of the initial set of N TLMs, and of course 



is the correct first threshold. Eq. (36) is just a scalar 



equation for the first mode intensity, yielding 



h 



1 



(i) 



(40) 



where Xii = Vi plays the role of the mode volume, en- 
hancing the power slope if mode one is more evenly dis- 
tributed over the gain volume. We should point out that 
I\ should be thought of as the intensity within the cav- 
ity. The emitted power is found by integrating the pho- 
ton flux associated with the TLM it M (r) over a surface at 
infinity |16j ; the transmissivity of the cavity is implicitly 
contained in the calculation of the TLM. In Appendix 
A we show that this power output can be related to a 
volume integral of the TLM over the gain region of the 
cavity, and that for single-mode lasing within the SPA- 
SALT one finds 
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= fei Jd d rF(r)\ Ul \ 2 
1 2tt Jci d rF(r)u2|ui| 2 



(A> 



(i) 



(41) 



Recently this equation was found to agree very well with 
the output power of a novel surface-emitting photonic 
crystal laser calculated using non-linear FDTD methods 

Using the procedure described above we now determine 
the mode v with the lowest interacting threshold and the 

(2) 

correct 2x2 matrix above this threshold. The second 
threshold is found to be 

,(2) Xll - X21 



-^O.int — 



Xll 



Dl 



iX21 



1-X 



-D, 



(2) 



where the interaction coefficient 



A 2 = 



D, 



(i) 



(2) 



X21 



> 0. 



Xii - X21 



(42) 



(43) 



Note that, as D}f' > as long as the modal in- 

teraction coefficient X21 is non-vanishing, the interact- 
ing second threshold is higher than the non-interacting 
threshold. The gain clamping limit is reached when 
A2 — > 1 => X21 — > XiiDq/Dq, and the first mode sup- 
presses any second mode for all values of the pump. One 
sees that strong overlap X21 ~ X11 leads to gain clamping 
as we expect. Also if the second mode has significantly 
lower Q-value or is away from the center of the gain curve, 
the ratio Dq /D^ is reduced leading to gain clamping 
for smaller values of X2i- One way to achieve this limit is 
in a microcavity laser with passive cavity modes spaced 
more widely than the gain bandwidth. 

(2) 

When the pump exceeds the second threshold D / nt , 
the modal intensities I\ and I2 are obtained from 



Eq. (36), 



h 



h = 



X22/D ( a 1) -Xl2/Df 

Ti(xiiX22 - X12X21! 



(D - D\ 



Xii/D { 2) -X2i/D { o 

r 2 (xnX22 - X12X21! 



(i) 



is given by 



where the modified intercept D\ 



'(i) 




1*P = 



X22 — X12 



X22 



(2) 



-D 



O.inf 



(44) 



(45) 



(46) 



fiyXi2 



The change in intercept indicates that the first mode in- 
tensity has a negative kink at the second mode thresh- 



old (L> (1) < £>o!it)> 

as can also be seen directly from 
the slope of Ii, which is reduced from its value of 
l/(r 1 x 11 £>Q 1 ^) in the interval below the second thresh- 
old. This kink is always negative because the turning 
on of a second mode reduces the slope efficiency of the 
laser in the first mode, but vanishes when the interac- 
tion coefficient \\2 ~^ and the two lasing modes act 
independently. 
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FIG. 6: (Color online) Modal intensity versus pump strength 
in a ID slab resonator. The description of the resonator is 
given in the caption of Fig. [3] Open symbols show the numer- 
ical solutions of ( 28 1 and s olid lines are the results of single- 
pole approximation (Eqs. ( |40|44)45 1). The color scheme is: 
blue (Mode 1), red (Mode 2), black (total intensity in the 
two-mode regime). 



F. Tests of the SPA-SALT 

To test the results derived above, we first revisit the 
ID laser studied in Scction |HI B| Fig.[6]shows the growth 
of modal intensities with Dq. In the single- mode regime, 
the result given by ( 40 ) agrees very well with the numer- 



ical solution of (28), indicating that the single-pole ap- 
proximation is almost exact. Consequently, the second 



threshold D, 



(2) 




0.892 is also accurately predicted by 

-,(2) 



(42), which gives Dq = 0.899. In the two- mode regime 
we still find good agreement, but the SPA-SALT slightly 
overestimates the suppression of the second mode. Nev- 
ertheless, the total intensity is in good agreement with 
the full SALT solution. 

To demonstrate the accuracy of the SPA-SALT in 
cases where the mode density is high, we study a uni- 
formly pumped 2D disk laser of radius R — 1 and in- 



dex n = 3.3 + 10" 



The gain is assumed to center at 



Re[nk a R] = 66 with width 7^ = fc a /40. Now there exist 
high-Q whispering gallery modes, and we find that first 
two thresholds are very closely spaced (see inset; Fig. [7]) 
and are four orders of magnitude smaller than those in 
the ID example just treated. The SPA-SALT correctly 
captures the intensity crossover of the first two modes 
shortly after the second one turns on, and its prediction 
for the first three modes remains impressively accurate, 
even after the onset of the 7th mode. As we have seen 
in Fig. [5jb), the higher order mode(s) are less single- 
pole like compared to the lower order ones. Thus we ex- 
pect the SPA-SALT not to work as well for higher order 
modes; this can be seen from the noticeable differences 
in the 5th (black) and 7th (cyan) thresholds given by the 
SPA-SALT and the full SALT results. Nevertheless, the 
slopes of all the higher-order modal intensities are still 
largely correct. 
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FIG. 7: (Color online) (a) Modal intensity versus pump 
strength in a 2D disk laser of uniform index n = 3.3 + 10~ 5 i, 
uniformly pumped. Squares show the numerical solutions of 
( 28 1 and solid lines are the results of the single-pole approx- 
imation ( |37[ ). Inset: Zoomed view near the first two thresh- 
olds, (b) Modal gain versus pump strength for the first ten 
TLMs, calculated with the full SALT, indicating that mode 
ten will never turn on due to modal interactions, as predicted 
by the SPA-SALT. The dashed line indicates the fully sup- 
pressed tenth mode. The first two modes are too close to- 
gether to be distinguished in this plot. Modal gain is defined 
in terms of eigenvalues of the modified lasing map and a mode 
reaches threshold when the modal gain reaches unity 26 . 



As noted, Eq. (39) gives a criterion for a complete 
suppression of modes after a certain number of modes 
have turned on. Typically if a mode is completely 
suppressed this equation gives a negative (unphysical) 
result. This happens for the tenth TLM in the current 
example. Indeed, the full SALT calculation, using the 
modified threshold matrix [26], confirms the prediction 
that mode ten will never turn on (see Fig. [7jb)). 

The SPA is better satisfied the less open is the laser 
cavity. The random laser is a system in which there is 
no conventional cavity, only multiple scattering to slow 
escape. In the most challenging case of a weakly scatter- 
ing RL, it has no sharp linear resonances at all, only the 
presence of the gain medium allows strong preference for 
certain frequencies [25]. In Ref. [25] the modal intensi- 
ties for a 2D RL were found within the full SALT theory 
to be a nonlinear function of the pump, unlike all other 
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FIG. 8: (Color online) Modal intensity versus pump strength 
in a random laser. The cavity is a disk of radius R, uniformly 
pumped with k a — 60/i? and 7x = 1/7?. The scattering mean 
free path is I = -R/3, significantly shorter than the systems 
studied in Ref. [25]. The solid curves give the exact SALT 
solution, and the dashed curves the SPA-SALT solutions. 



cases studied. Thus we do not expect the SALT to apply 
there. Even when the disorder scattering is increased in 
the RL in order to increase the Q, and the intensities are 
linear in the pump, we find that the SPA-SALT, while 
it still gives good qualitative results, does not give good 
quantitative agreement with the exact SALT solutions, 
as shown in Fig. [8] 



V. SUMMARY AND CONCLUSIONS 

We have presented an improvement of steady-state ab 
initio laser theory based on the TCF basis that allows 
one to solve the self-consistent SALT equations more 
efficiently for resonators which are spatially inhomoge- 
neous (as is usually the case) and/or with inhomogeneous 
pumping. This completes the development of the ab- 
initio theory based on the stationary inversion approxi- 
mation, originally proposed in 2006 [16] and improved in 
several subsequent papers [17 ] I23H26] . This theory takes 
into account the openness of the cavity exactly in terms 
of TCF or UCF basis states, and includes the nonlinear 
hole-burning interactions to infinite order. Besides pre- 
dicting interacting thresholds and intensities, the theory 
captures subtle effects such as the change in shape of the 
lasing modes, and the variations in their frequencies as 
the pump is increased well above threshold. 

Using the TCF basis and the single-pole approxima- 
tion |16j , we have derived a simplified version of the the- 
ory, the SPA-SALT, which predicts a linear increase of 
all lasing intensities. The relevant slopes and interacting 
thresholds can be found with negligible computational 
effort, once the linear problem of the non-interacting 
threshold lasing modes is solved. Explicit analytic so- 
lutions were for the few-mode lasing regime, illustrating 
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important qualitative features of multimode lasing with 
modal interactions. In particular, an analytic condition 
was found for the "gain-clamping" transition, in which 
higher modes are completely suppressed by modal inter- 
actions. For non-trivial examples, the SPA-SALT agrees 
well with the full SALT, although its breakdown for very 
low-Q systems such as random lasers was also found. Al- 
though further work is needed to determine the regime of 
quantitative validity of the SPA-SALT, there is already 
evidence that it will be possible to dramatically improve 
the modeling of realistic and complex laser structures in 
two and three dimensions |39) , by reducing the nonlin- 
ear lasing computation to almost the same level of diffi- 
culty as the linear problem of finding the threshold lasing 
modes. Even when the SPA-SALT is not a good approx- 
imate theory, the full SALT equations in the TCF basis 
will improve steady-state lasing calculations by many or- 
ders of magnitude compared to brute force time-domain 
simulations. 
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convenience: 

-p = 

Z7T 



d J r 



l + h(r) 



Im[ e (r)] H*„(f)| s 



(A4) 

This result states that the total power radiated by each 
lasing mode equals the power that the gain medium de- 
livers into that mode, minus the power that the mode 
loses through material dissipation (described by Im[e]). 

It is instructive to consider the modal power in the 
single-pole approximation. Let us suppose that Im[e] = 
0. Combining the general expression for in (A3) with 
the SPA-SALT expression * M ; 



I„ Uu, we obtain 



M ~~ 47Tfc, 



72 < 



c.c. 



d 2 r [u M V 
|r M £> %|d 2 rF(r)| V | 



(A5) 
(A6) 



Here we have used ( 19 ) to express Im^] in terms of the 
SPA laser thresholclT>Q . As noted in the main text, in 
the single-mode regime (p — I), the modal power has a 
particularly simple form: using (33), we can write 7 M in 
terms of the pump Dq, to obtain 



Pi = 



k x Jd 2 rF(r)\ Ul 



2tt $d 2 rF(r)ul\ Ul \ 2 



{Do -D?) 



(A7) 



Appendix B: Comparison to Mandel approach 



Appendix A: Modal output power 

In this appendix we derive the modal output power 
of a 2D cavity from the (internal) modal intensity. The 
output power is the total flux of the Poynting vector, 
taken across a loop T enclosing the cavity: 



P = 



1 

-in 



ds i 



Ex B 



(Al) 



In Q , the (out-of-plane) electric field is written as a sum 
over the modal fields \&„(r), and a similar expression may 
be written for the (in-plane) magnetic field. We find the 
time-averaged total output power (V) = Y^u,Ppi with 
the modal power P^ given by 



47rfc 



- (fdsh- [^Vtf* -c.c] (A2) 



ink 



fi Jc 



d 2 r [«VV 2 ^ - c.c] . (A3) 



In the last step we have used the Gauss' law. Here ^ 



and 



and P^ are measured in their natural units e L 
introduced when deriving Eq. Q. Using Eq. (A3), to- 
gether with the wave equation |9| and its complex con- 
jugate, gives us Eq. (11), which we reproduce here for 



In Refs. [22 [53] Mandel and coworkers treated the 
infinite-order modal interactions in a Fabry-Perot cav- 
ity in the single-mode and two-mode regimes. They used 
the approximations of stationary inversion, and pump- 
independent lasing modes and frequencies, similar to 
the SPA-SALT (the SALT of course includes the pump- 
dependence of the lasing modes and frequencies [2"5]). 
Unlike the SPA-SALT, they assumed that the fixed lasing 
modes were hermitian closed cavity modes (sine waves of 
real wavevector). They did not derive a version of the 
basic constrained linear equation p3| of the SPA-SALT, 
but instead they derived a single-pole closed-cavity ver- 
sion of Eq. (28). For the single-mode case, Eq. (5) of 
Ref. [25] is of exactly the same form as Eqs. (43), (54) 
of Ref. [TH], the earliest version of the SALT, except for 
their use of closed cavity modes. Ref. [H] applies the 
single-pole approximation to the direct map but treats 
the openness of the cavity exactly using non-hermitian 
constant flux states; this approximation is not exactly 
equivalent to the SPA-SALT, which uses the SPA on the 
inverse map, but gives very similar results to the SPA- 
SALT at large pump strength. 

It is interesting to compare the two methods for the 
simple case of a uniformly pumped ID dielectric slab 
laser of the type considered in [16l [TT], [24] (see inset, 
Fig.0. We will compare Mandel's approach to the full 
SALT, the most complete form of our theory. Thus our 
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approach differs from Mandel in two major ways. First 
we take into account the openness of the cavity exactly 
and second we allow for the change in the lasing modes 
and modal frequencies above threshold. To vary the qual- 
ity factor of the cavity, we choose four sets of parame- 
ters {n, k a } = {1.5, 40}, {3, 20}, {5, 20}, and {10, 20} (k a 
is the frequency of the gain center). We have shown 
in Ref. [T7] that for the first two sets of parameters the 
SALT and numerical solutions of the MB equations agree 
very well, so we can take the SALT results as correct. 

The rescaled model intensity (1'^ = L^/^) in the single- 
mode regime in Mandel's approach is given in our nota- 
tion by 



0.15 



>- 0.1 - 




I' (Do) 




(Bl) 



is produced using Eq. (Bl I, which has no dependency on the 



FIG. 9: (Color online) Rescaled intensity 1'^ = r M J M of the 
first mode in a ID slab resonator. The cavity is open on both 
sides and the pump is taken to be spatially uniform. Do is the 
The dependence on the refractive index of the cavity is pump intensity and. DqJ_ is its threshold value. The solid line 

contained in the first threshold, which is not calcu- 

lated in the Mandel approach, but is assumed known and 
used to normalize the pump. The gain parameters (k a 
and 7jJ only enter in the scale factor (r^) and implicitly 

again through D^p . Note that the Mandel single-mode 
result has an additional square root dependence on the 
pump, which is not present in the SPA-SALT. This dif- 
ference arises because, as already noted, the single-pole 
approximation is made at a different point in the two 
derivations. The full SALT theory does not predict a 
universal linear dependence on pump and indeed for very 
low-Q lasers, such as random lasers, the dependence can 
be non-monotonic |25j . 

In Fig. [9] we compare the result given by Eq. ( Bl I to the 
SALT. As one might have expected, the two approaches 
agree well for the higher Q cases (n = 10, 5) but a sig- 
nificant disagreement in the slope of the intensity curves 
appears for the lower Q (n = 1.5,3) cases. Nonethe- 
less, the Mandel approach for the single-mode case is 
qualitatively better than HS, which shows an unphysical 
saturation [I"6l IT7]. 

Next we compare the value of the interacting second 

(2) 

threshold Dq ( nt given implicitly in Mandel's method by 



cavity index or length. The four dashed lines are the results of 
the SALT with different cavity indices and atomic transition 
frequecies. The upper bound of single-mode lasing in the 
n — 10 (high-Q) case is near Do/Dq 1 ~ 1.15, and the intensity 
overlaps with the solid curve. As the cavity index/Q-factor is 
reduced, the black curve differs more and more from the SALT 
result, whose accuracy has been proven by comparing with the 
time-dependent simulation of the steady-state solutions of the 
MB equations (Ref. [17]). 



qualitatively better than that of HS but is not as accu- 
rate as the SALT and the SPA-SALT, both of which are 
based on general computational algorithms applicable to 
arbitrary cavities. 



2.2 



-° 1 



D, 



(i) 



D 



(2) 



2) 



D, 



(2) 



4^ 



« (D. 



(2) 



iio 
Q 
II 



D, 



(2) 



D, 



(2) 



(B2) 

and the result of the SALT in the four cases listed above. 
We find that Mandel's approach consistently underesti- 
mates the strength of the modal interactions and deviates 
relatively little from the non-interacting threshold values 
(see Fig. [To] ) . The highest Q case agrees most closely 
with the SALT, but there is some non-monotonic behav- 
ior of the thresholds with Q value in the SALT which we 
did not analyze in detail. We conclude that the effect 
of openness accounts for the main difference between the 
SALT and the Mandel approach, in a Fabry-Perot cav- 
ity in which both can be applied. Mandel's approach is 
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FIG. 10: Second interacting threshold in a ID slab resonator. 
The cavity and parameters used are the same as in Fig. [9] 



,(-') 
'0 

interaction 



Dq is the second threshold value in the absence of modal 
The solid line and crosses are the solution of 



Eq. ( B2 1 , and the dotted line indicates the non-interacting 
case (7 = 1/w). The results of the SALT are indicated by the 
different symbols explained in the legend. 
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Appendix C: Perturbative calculation of corrections 
to the SPA-SALT 



The major approximation in the SPA-SALT is re- 
placing the expansion (26) with a single term, = 



a^Ufj,. In this appendix we derive the first order 
expression for the non-dominant expansion coefficients 
in the single-mode regime. Assuming the dominant 



component is a 
h\ Ul (r)\ 2 . Eq. i 



we approximate h(r) by |aiMi(r)| 2 
is then 



30 1 for a. 



FIG. 11: Expansion coefficients of the first lasing mode (left) 
and the second lasing mode (right) at the second threshold 
in a ID slab resonator. The solid curve is the solution of 
( |28[ ) and the dashed line is given by the approximation ( C2 1 
and 



( |C3[ ), respectively. Notice that the vertical axis is in the 
logarithmic scale; the expansion of the first/second mode is 
dominated by the first /second UCF state with a weight of 
90%/84%. 



n(n> 1) 



M> 0-n — T HI \1\ — a m, 



(CI) 



A » rn A 



By inserting the expression ( 40 1 for Jj , derived in the 
single-pole approximation, into Eq. (CI), we reduce the 
latter to a set of inhomogeneous linear equations of 
a n(n>i)- Eq. (CI) can be further simplified by keeping 
only the a\ term m the sum, which leads to 



We can also derive an analytical expression to evaluate 
the non-dominant expansion coefficients of the second 
mode when it turns on. We assume that its dominant 
component is 0,% ~ , and derive 



Xnll 
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Note that 



u ' 



which is 
to 



^n( fc i )/7i is not 
?7n=^Wi/7Ai- In Fig. I Ilia) we compare (C2) 
the numerical solution of ( |28[ ), and they agree very 
well even in the logarithmic scale. The system is the 
inhomogenous ID resonator considered in the main 
text, and the pump is tuned to the second threshold 
(D = 0.892). 
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72^0 



in the same way (C2) is derived, 
using (36) that the ratio becomes 1 when 
should 



It is easy to check 
n = 2 as it 



The result above is compared with the multi- 
pole expansion ( 28 ) in Fig. fTTVb) . 
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